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Abstract 
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. We present numerical studies of the transport dynamics of a charged soliton in conjugated 

O ^ 

\ polymers under the influence of an external time-dependent electric field. All relevant electron- 

> 

Q ' phonon and electron-electron interactions are nearly fully taken into account by simulating the 

monomer displacements with classical molecular dynamics (MD) and evolving the wavefunction 
, for the vr electrons by virtue of the adaptive time-dependent density matrix renormalization group 

^ ^ (TDDMRG) simultaneously and nonadiabatically. It is found that after a smooth turn-on of the 

■ external electric field the charged soliton is accelerated at first up to a stationary constant velocity 



as one entity consisting of both the charge and the lattice deformation. An ohmic region (6 mV/A 
< Eq < 12 mV/A) where the stationary velocity increases linearly with the electric field strength 



' is observed. The relationship between electron-electron interactions and charged soliton transport 

^ ■ is also investigated in detail. We find that the dependence of the stationary velocity of a charged 

soliton on the on-site Coulomb interactions U and the nearest-neighbor interactions V is due to 

(N 

' the extent of delocalization of the charged soliton defect. 
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I. INTRODUCTION 



For more than three decades, there has been a great number of both academic and 
industrial interest in irans-polyacetylene (PA), since the discovery that its electrical con- 

nnn 

ductivity can be improved greatly through charge doping or photo-excitation pj, 12, 13[. It 
is now widely accepted that the enhanced conductivity of trans-PA is due to its inherent 
tendency to show solitonic behavior. Because there exists a twofold degeneracy of ground 
state energy in trans-PA distinguished by the positions of the double and single bonds, the 
soliton excitations take the form of domain walls separating different districts of opposite 
single and double bonds alternation patterns. The reduction or the oxidation of the trans- 
PA chain will result in polaron defect and/or spinless charged soliton defect. At least at 
low doping levels, conductivity of tmns-PA is thought to result from the motion of these 
quasiparticles. [J, l5|, la, 

mm 

Obviously, it is of great importance to understand the underlying mechanisms of charge 
carrier transport under an external electric field in order to modulate or devise new functional 



materials based on trans-PA. There have been extensive theoretical studies o 
of solitons in trans-PK under the influence of external forces 
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the dynamics 



6J] based on the 



Su-Schrieffer-Heeger (SSH) model , an improved Hiickel molecular orbital model, in 
which TT electrons are coupled to distortions in the polymer backbone by the electron-phonon 
interaction. In order to improve the situation that there are no electron-electron interactions 
included in SSH calculations, Forner et al. [l^ performed Pariser-Parr-Pople (PPP) 
simulations for the soliton dynamics in trans-PK at the self-consistent field (SCF) level and 
found that the soliton width decreased with the inclusion of electron-electron interactions 
and that the soliton velocity {vs) is smaller than that found in SSH calculations. By using an 
ah initio molecular dynamic with a single spin-unrestricted determinant. Champagne et al. 
characterized the motion of carbon and hydrogen atoms in a small polyacetylene chain with 
a positively charged soliton defect driven by an external electric field and found that the 
motion consists of dominant longitudinal CH displacements combined with less prominent 
H wags and CH stretches. 151] 

However, recently a lot of theoretical calculations of the static properties of conju- 
gated polymers have shown that the electron correlation effect plays a very important 
role in determining the behavior of the charge carriers in conjugated polymers, [si Q, 16, 
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29( 1 For example, it is found that the soh- 



ton width depends strongly upon the inclusion of the electron correlation effects and the 
electron correlation effects are substantial for characterizing the geometrical and charge 
features of vr-conjugated chains bearing a charged soliton defect without counterion. 28] 
Therefore, theoretical simulations for soliton dynamics with electron correlations consid- 
ered, which step beyond the SCF level, are highly desirable. However, the calculations 
for large polymer systems with traditional advanced electron-correlation methods such as 
configuration interaction method (CI), mult i- configuration self-consistent field method (MC- 
SCF), many-body MoUer-Plesset perturbation theory (MPn) and coupled cluster method 
(CC) are currently still not feasible due to the huge computational costs. Fortunately, 
the density-matrix renormalization group (DMRG) method, firstly proposed by White 



in 1992 



30, 



31 



32l |. can be used instead: it has been demonstrated and widely ac- 



cepted that at economic computational cost the accuracy of DMRG is very close to ex- 
act full configuration interaction (FCI) calculations for one-dimensional (ID) or quasi- 
ID systems, where DMRG is the most performing numerical method currently available. 
Since Pang and Liang introduced the DMRG method into the studies of conjugated poly- 
mers with an alternating Hubbard model in 1995 Q], a lot of interesting papers by 
many groups focusing on the static investigations of conjugated polymers with the ap- 



)lication of DMRG within various einpirical or semi-ernpirical models have 
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In this paper, combining the SSH and the extended Hubbard model (EHM), we simulate 
the motion of the charged soliton in trans-PA under an applied external electric field based 
on a nonadiabatic approach, in which we evolve the vr-electron wavefunction using a time- 
evolution scheme for DMRG (the adaptive time- dependent density matrix renormalization 



group (TDDMRG) 58|, |59|) and simulate the motion of the nuclei in the lattice backbone 



by classical molecular dynamics (MD) under the couplings with the vr-electron part simul- 
taneously, denoted for convenience as TDDMRG/MD. The aim of this paper is to give an 
exhaustive picture of charged soliton transport in trans-PA at a theoretical level with all 
relevant electron-electron interactions and correlations included and to show how the on- 
site Coulomb interactions and nearest-neighbor electron-electron interactions influence the 
behavior of charged soliton transport in trans-PA. 
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II. MODEL AND METHODOLOGY 



A. Model 

We use the well-known and widely used SSH Hamiltonian combined with the ex- 

tended Hubbard model (EHM), and include the external electric field by an additional term: 

H{t) = Hel + HE(t) + Hiatt (1) 

This Hamiltonian is time-dependent, because the electric field E{t) is controlled explicitly 
by time. 

The TT-electron part includes both the electron-phonon and the electron-electron interac- 
tions, 

Hel = "~ ^n,n.+l(c^+l,crCn,o- + h.C.) 

n,a 

+ ^(c^.CT^n.CT - -;^{ct-a(^n~a " 2) (2) 



where tn,n+i is the hopping integral between the n-th site and the (n + l)-th site, while 
U is the on-site Coulomb interaction and V denotes the nearest-neighbor electron-electron 
interaction. Because the distortions of the lattice backbone are always within a certain 



limited extent, one can adopt a linear relationship 



3etween the hopping integral and the 



lattice displacements as t„,n+i = — — Un) sl, O], where to is the hopping integral 

for zero displacement, m„ the lattice displacement of the nth CH monomer, and a is the 
electron-phonon coupling. 

Because the atoms move much slower than the electrons, we treat lattice backbone clas- 
sically with the Hamiltonian 

Hlatt = Y '^i^n+l - Unf + y J] ^« ' 

n n 

where K is the elastic constant and M is the mass of a CH monomer. 

The electric field E{t) directed along the backbone chain is uniform over the entire system. 
The field which is constant after a smooth turn-on is chosen to be 

E,exp[-{t-Tcf/T^l for t<Tc, 

(4) 

for t > Tc, 




where Tc, T^r and Eq are the center, width and strength of the half Gaussian pulse. This 
field gives the following contribution to the Hamiltonian: 



^Eit) = \e\ ^{na + Mn)(c+^c„,a - ^)^(^) 



(5) 



where e is the electron charge and a is the lattice constant. The model parameters are 
those generally chosen according to Su et aVs pioneer work |9|: to=2.5 eV, a=4.1 eV/A, 
K=21 eV/A^, M=1349.14 eVfsVA^ a=1.22 A. We notice that there are some debates on 
the accuracy of some of the parameters. For example, Ehrenfreund, et al. proposed to 
use K=AQ eV/A^ instead of 21 eV/A^ to fit the resonant Raman scattering experiment 
60|, and an additional linear term was suggested to be included into to Eq. [3] for the 
more proper description of the elastic energy of the C-C a backbone oj] since a is assumed 
to be the equilibrium lattice spacing of the undimized chain, including both a- and vr- 



3onding. 
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ur the purpose of making contact to other theoretical studies on this problem 
"691 . we choose the current parameter schemes. In order to 
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prevent the contraction of the chain, which might result from our simplified lattice elastic 
energy form, we supplement a constraint of fixed total length of the chain. The values of Tc 
and Tw (Tc=30 fs and Tvi/=25 fs) are taken from the paper of Fu et al^\ 



B. the TDDMRG/MD method 



Before moving to the dynamical simulation, we determine the static electronic structure 
of the system in the absence of the external electric field by virtue of DMRG. The DMRG 



method 



30 



3l| has emerged over the last decade as the most powerful method for the 



simulation of stro ngly correlated ID quantum systems. We would like to refer readers to 



recent reviews 



32 



Id 



71| for the technical 



we utilize energy minimization methods 72 



detai 



s of DMRG. During the DMRG calculations. 



73| to optimize the initial lattice configuration 



{un\ for the chain backbone with the energy gradients being calculated according to the 
Hellmann-Feynman theorem 0, [t^ . 

After the lattice configuration and the vr-electron wavefunction for the initial state are set, 
we turn on the external electric field smoothly and started the real-time simulation of the 
dynamic process of the charged soliton transport along the model chain by the nonadiabatic 
TDDMRG/MD method. The main idea of the TDDMRG/MD method is to iteratively 





Step 0: 

Perform energy minimization of Hff"") by DMRG, 

get atoms initial positions 1^*="' and initial electron wavefunction ipc^"' , 

choose short At 








Step 1: 

Perform time simulation for electron evolution by TdDMRG, 
get new electron wavefunction ijjC+'i') 










Step2: 

Get forces F=-VU(r) and a= F/m 








Step 3: 

Move atoms: r('+*') = i<')+ vW At+ a At V2 








Step 4: 

Move time forward: t=t+At 





FIG. 1: Procedures for the TDDMRG/MD method. 



repeat the procedures of evolving the yr-electron part by the adaptive TDDMRG and moving 
the backbone part by classical MD simultaneously and nonadiabatically. The flowchart of 
the procedures for the TDDMRG/MD method is illustrated in Fig. [TJ 



The adaptive TDDMRG method 



58 



591] , which is an extension of standard DMRG using 



Vidal's time-evolving block-decimation (TEBD) algorithm [76|, is a recently developed very 
efficient numerical method to solve the time-dependent Schrodinger equation 

^{t + At) = U{At)-^{t) with U{At) = exp{-iH{t)At/h) . 

Time evolution in the adaptive TDDMRG is generated using the Trotter-Suzuki decomposi- 
tion of the time-evolution operator U (At) of Eq. ([U]). Since the Hamiltonian operator for the 
TT-electron part of Eq. can be decomposed into a sum of local terms hn that live only on 
neighboring sites n and n + 1, U{At) can be approximated by an nth-order Trotter-Suzuki 
decomposition [t^], e.g., to second order, 

UiAt)= H f/„(At/2) H Un{At) W U4At/2) + 0{At') (7) 

n£even nGodd nGeven 

The Un{At) are the infinitesimal time-evolution operators exp{—ihn{t)At/h) on the bond 
linking sites n and n + 1. The ordering within the even or odd products doesn't matter, 
because "even" and "odd" operators commute among themselves. Then we can incorporate 
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the applications of these operators t/„(At) successively to some state into finite-system 



DMRG sweeps 



3Q|, 



31 



32| . Each operator Un{At) is applied at a finite-system DMRG step 



with sites n and n + 1 being the active sites, i.e., where sites n and n + 1 are represented 
without truncation. For example, in a L-site ID chain, at a finite-system DMRG step with 
n and n + 1 being the active sites, we can use Eg. ([8]) to describe a DMRG state in the form 
of matrix product states (MPS) 78, 79, 80, 8l|, 82| (also known as finitely correlated states), 

|^)= A'[a,]A'[a,]...A'^[aL]\a^...aL), (8) 

(Tl...(7i 

where denotes the state space for site i, matrices A* are for site i, and the first and 

last A-matrices are taken to be row and column vectors. Then it is easy to describe the new 
state after the operation of f/„(At) 

(9) 

without any additional error, because f/„(At) only acts on the part of Hilbert space 
iWn)Wn+i)) which is exactly represented. Similar to conventional DMRG, in order to prevent 
the exponential growth of matrix dimension for continuing the finite-system sweep, DMRG 
truncations must be carried out. But differently, the adaptive TDDMRG uses ?7„(At)|\i/) as 
the target state instead of |^) to build the reduced density matrix. At the end of several 
sweeps, when all the local time-evolution operators have been applied successively to "^(t), 
one can get the vr-electron wavefunction \E'(t + At) for the new time t + At, and continue to 
longer times. In the context of ID electronic and bosonic systems, the adaptive TDDMRG 
has been found to be a highly reliable method, for example in the context of magnetiza- 



tion dynamics 83|] of spin-charge separation 8J, l85|] or far-from equilibrium dynamics of 



ultracold bosonic atoms 



86|. 



During the TDDMRG/MD simulations, after we evolve the vr-electron wavefunction from 
ttot + At and get \E'(t + At) using the adaptive TDDMRG technique, we start classical MD 
simulations for the motion of atom nuclei in the chain backbone. The force acting on the 
nth site is equal to the negative energy gradient of nth site's displacement. Here, we take 
the non-Born-Oppenheimer effect into account and calculate these forces in a nonadiabatic 
way, as in Eq. ffTOj) . 

Mun{t) = -^ = -{m\^\m) (10) 



dUr. 



dUr 



In the case of our SSH-EHM Hamiltonian, we can calculate these forces as 

MUn{t) = - K[2Un - Un+l - ^n-l] 

+ \e\E{t)[l - (M/(t)|5^<,c„,.|vl/(t))] 

CT 

+ a5^[(^(t)|c++i,,c„,^ + /i.c.|^(t)) - {^>{t)\cl^Cn^i,. + h.c.\^{t))]. 

(7 

Once we obtain the forces acting on the atom nuclei, we can simulate the movement of CH 
monomers according to Newton's laws of motion. 

When the new lattice configuration {m^} for the backbone chain is determined, we move 
the time forward and go back to the TDDMRG step to evolve the vr-electron wavefunction 
and then repeat this process as long as we need or want. Upon repeating this process again 
and again, we can simulate a long time vibration process of the chain backbone, while the 
TT-electron part is also well simulated quatum mechanically in real time by virtue of the 
adaptive TDDMRG. 



III. RESULTS AND DISCUSSION 



To investigate the dynamic properties of solitons, we simulate the charged soliton trans- 
port process in a single model conjugated polymer chain under a uniform external electric 
field. In order to present the twofold degeneracy of the ground state energy in soliton 
defect, the finite-sized model chain is specified to have an odd number of sites. In our calcu- 
lations, we simulate a positively charged model chain containing = 101 CH monomers and 

= 100 TT electrons. Due to the inherent electron-hole symmetry, our employed parame- 
terized SSH-EHM model can not be used to distinguish positively and negatively charged 
defects. In order to simulate a longer dynamic process of charged soliton transport within 
the chain with a certain length, the charged soliton soliton is initially located in the left side 
of the chain instead of being centered in the middle of the chain as usual. The charged soli- 
ton is initially centered around site 21 through imposing a constraint of reflection symmetry 
around site 21. After then, a 160-femtoseconds (fs) dynamical process, which is generally 
long enough for the charged soliton to be transported from the left end to the right end in 
our model chain, is simulated by the nonadiabatic TDDMRG/MD method. 

Because this is the first time that we perform the dynamical simulations using MD com- 
bined with TDDMRG, we will firstly give a demonstration that one can achieve reliable 
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FIG. 2: The net charge of a charged sohton as a function of the site index at t = 160 fs calculated 
with different values. (At=0.05 fs, ^o=6.0 mV/A, U=2.0 eV, V=1.0 eV) 



results for ID correlated quantum systems with appropriate values chosen for two key pa- 
rameters. After the submission of this work, we notice that very recently Zhao et al have 
applied such a TDDMRG/MD method to the studies of polaronic transport within a com- 
bined model of SSH and Hubbard model and the excellent agreement of their TDDMRG/MD 
results with the exact numerical results for a noninteracting chain (f/=0, V=Q) verified the 
reliability of the TDDMRG/MD method. 68|| Therefore, here we just mainly focus on show- 
ing the convergence behavior of the numerical results with different values for the two key 
parameters and showing how to choose suitable parameter values for the purpose of achiev- 
ing reliable results. The two most important parameters in the nonadiabatic TDDMRG/MD 
simulations are the DMRG truncation weight and the time step At. 

Fig. [2] shows the net charge distribution pictures at t = 160 fs calculated by different 
DMRG truncation weight values. As we have mentioned in the above section, is a 
key parameter for estimating the DMRG truncation error of the wavefunction and local 
quantities. So, should be also very vital in our dynamical simulations by MD combined 
with the adaptive TDDMRG. As can clearly seen from Fig. [2], a big ep value will lead to large 
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FIG. 3: The net charge of a charged sohton as a function of the site index at t = 160 fs calculated 
with different At values, (e^ = 5.0 x 10-^ -£0=6.0 mV/A, U=2.0 eV, V=1.0 eV) 

deviations from the results obtained by small ep values, and even the peak of the net charge 
density may be shifted to a totally wrong position. Meanwhile, when we decrease value, 
the net charge values converge gradually; when ep has been small enough (ep < 1.0 x 10~®), 
the net charge values will not improve much if we decrease ep value further. This situation 
of how ep values influence the final simulated results is very similar to that in standard static 
DMRG calculations and it implies that our dynamical simulations may yield reliable results 
close to the real case if we control to be small enough (ep < 1.0 x 10~^). 

We show the net charge distribution pictures at t = 160 fs calculated by different At 
values in Fig. [3l It is clearly shown that because of large Trotter error the too large At 
value (0.0625 fs) results in apparently different net charge distribution pictures comparing 
to the other At values, while the At values of 0.05 fs, 0.02 fs and 0.01 fs produce similar 
charge distribution patterns which are at least qualitatively in good agreement with each 
other although there are still some slight differences. Different from the traditional time 
evolution methods, a smaller time step At doesn't always lead to higher accuracy in the 
adaptive TDDMRG simulations. Previous studies have elucidated that there are two main 
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FIG. 4: Time evolution of net charge distribution in a trans-PA chain with the transport of a 
charged soUton under an external electric field. (The charged soliton moves from the left to the 
right.) (At=0.05 fs, Cp = 5.0 x 10"^ ^o=6.0 mV/A, U=2.0 eV, V=1.0 eV) 

error sources in the adaptive TDDMRG: the Trotter error and the accumulated DMRG 



truncation error. 



83 



87| Although a smaller time step At will of course result in a much 



smaller Trotter error, it requires more DMRG iterations and accordingly it will lead to a 
larger accumulated DMRG truncation error. So, one should choose the time step At values 
carefully with a well-balanced account of these two sorts of errors. Therefore, we consider 
that a time step At value around 0.01~0.05 fs is reasonable for our current cases. 

We now move on to give the general time evolution picture of charged soliton transport 
in trans-PA. The time evolution of net charge distribution and the staggered bond order 
parameter r„ = (— 1)"(2m„ — m„+i — m„,„i)/4 in the singly doped trans-PA chain are shown 
in Fig. m and Fig. It can be clearly seen that, both the net charge and the geometrical 
distortions don't localize at one CH monomer. On the contrary, the defect is dispersed to 
several monomers around the soliton center. Actually, defect delocalization is an inherent 
feature of a soliton and very important for the solitonic transport. One can also clearly 
see that during the entire time evolution process the geometrical distortion curve and net 
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FIG. 5: Time evolution of the staggered bond order parameter r„ in a trans-PA chain with the 
transport of a charged soHton under an external electric field. (The charged soliton moves from 
the left to the right.) (At=0.05 fs, €p = 5.0 x 10"^ Eo=6.0 mV/A, U=2.0 eV, V=1.0 eV) 

charge distribution shape for the charged sohton defect always stay coupled and show no 
dispersion. This implies that the soliton defect is an inherent feature and the fundamental 
charge carrier in trans-PA. In Fig. |5l we can also see that a long lasting oscillatory "tail" 
appears behind the soliton defect center. This "tail" is generated by the inertia of those 
monomers to fulfill energy and momentum conservation; they absorb the additional energy, 
preventing the further increase of the polaron velocity after a stationary value is reached. 
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From Fig. [Hand Fig. [5], one can also find that, after the external electric field is turned 
on, the charged soliton is accelerated at first and then moves with a constant velocity as 
one entity consisting of both the charge and the lattice deformation. Fig. [U] of the temporal 
evolution of soliton center position and Fig. [7] of the temporal evolution of soliton velocity 
clearly support this observation. Therefore, only the stationary velocity will be considered 
for charged soliton transport in the following discussions. 

The dependence of the stationary velocity of a charged soliton Vg on the electric field 
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FIG. 6: Temporal evolution of the center position of a charged soliton under an external electric 
field. (Ai=0.05 fs, Ep = 5.0 x 10"^ Eq=Q.{) mV/A, U=2.0 eV, V=1.0 eV) 
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FIG. 7: Temporal evolution of the velocity of a charged soliton under an external electric field. 
(At=0.05 fs, €p = 5.0 X 10-^ Eq=6.0 mV/A, U=2.0 eV, V=1.0 eV) 

strength is shown in Fig. [HI We find that Vg increases with increasing electric field strength. 
Moreover, a very interesting behavior of Vg as a function of the electric field strength is 
observed. An ohmic region where Vg increases linearly with the field strength is found. This 
ohmic region extends approximately from 6 mV/A to 12 mV/A for the case of [7=2.0 eV 
and 1^=1.0 eV. Below 6 mV/A the soliton velocity increases nonlinear ly. The saturation of 
Vg with increasing the electric field strength is not found. 

In order to study the influence of electron-electron interactions on soliton transport, we 
focus on the stationary velocity of a charged soliton Vg calculated with different electron- 
electron interactions under the constraint that the other parameters are fixed. As should 
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FIG. 8: The stationary velocity of a charged sohton as a function of the electric field strength. 
(Ai=0.05 fs, ep = 5.0 x 10"^ U=2.Q eV, ^=1.0 eV) 
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FIG. 9: The stationary velocity of a charged soliton as a function of different U values. (At=0.05 
fs, Ep = 5.0 X 10"^ ^0=6.0 mV/A, V=0.0 eV) 

be noticed that, the real conjugated polymer is with weak couplings. Strong couplings will 
lead to too strong charge polarizations which are unrealistic. Considering that we are only 
focusing on the study of real conjugated polymer system, in this work we adopt only the 
weak-coupling parameters {U <2t,V <t). 

Firstly, we study the Hubbard model with only on-site Coulomb interactions f/, i.e., 
V = Q in the EHM. The dependence of the soliton stationary velocity Vg on the on-site 
Coulomb interactions U is displayed in Fig. [9l Interestingly, Vs is non-monotonic in f/: it 
grows to a shallow maximum at ?7 ~ 0.9 eV. Actually, the increase or decrease of Vg is 
strongly related to the degree of delocalization of the soliton defect. In Fig. [10] we show the 
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FIG. 10: The net charge of a static charged sohton as a function of the site index calculated with 
different U values, (e^ = 5.0 x 10"^ V=0.0 eV) 



net charge distribution pictures of a static charged soliton calculated with different U vaslues. 
Clearly, there are large charge fluctuations within the defect central part due to the existence 
of a soliton. While the on-site repulsion U increases, a large charge fluctuation becomes 
energetically unfavorable and the net charge has the tendency to be equally distributed at 
different CH monomers. Therefore, it can be clearly seen from Fig. [TU] that, the charge 
oscillation is reduced substantially with the increase of U. In fact, at the early stage of 
the increase of the on-site Coulomb interaction U (when the values of U are not too large), 
this increase is in favor of the transport of the soliton because the on-site repulsion has 
the tendency to push the electron (or hole) in the accumulated large charge density of the 
central part of a charged soliton to hop more easily to the neighbor site. In other words, 
within this stage the delocalization of the soliton defect is strengthened and consequently 
soliton transport becomes easier. Therefore, the stationary velocity of a charged soliton Vg 
increases while the value of U increases from 0.0 eV to 0.9 eV as displayed in Fig. [9l On 
the other hand, when the on-site Coulomb interactions U increase further and become to 
be dominant, the lattice tends to be occupied by one electron per site. This means the 
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FIG. 11: The staggered bond order parameter of a static charged soUton for several different U 
vaslues. (cp = 5.0 x 10"^ V=0.0 eV) 

electrons are more localized and accordingly the free motion of the charged soliton will be 
greatly restricted by large on-site repulsion U , so the delocalization of the soliton defect is 
decreased and consequently the transport of the soliton becomes more difficult. Therefore, Vs 
decreases while the value of U increases further from 0.9 eV to 4.0 eV as shown in Fig.[9l The 
change of delocalization level of the charged soliton defect can be directly viewed through 
the geometrical picture of the static soliton defect. In Fig. [TT], the staggered bond order 
parameter r„ of a static charged soliton calculated by different U values is shown. The 
geometrical distortion of soliton defect from the regular single and double bond alternation 
is normally described by a hyperbolic tangent function,!^, 9] 

'n-C' 



tanh 



L 



(12) 



where, Too is the staggered bond order parameter in infinite trans-PA chain without defect, 
C denoting the defect center, and L is the half-width of the soliton, determined in the way 
that it will give the best fitting of r„ against n. We extract the soliton width L as: 7.95 
{U = 0.0 eV), 8.18 {U = 0.9 eV) and 6.64 {U = 4.0 eV). These data certainly show that the 
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FIG. 12: The stationary velocity of a charged soliton as a function of different V values. (At=0.05 
fs, Ep = 5.0 X 10-^ £^0=6.0 mV/A, U=0.0 eV) 
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FIG. 13: The net charge of a static charged soliton as a function of the site index calculated with 
different V values, (cp = 5.0 x 10"^ U=0.0 eV) 

delocalization of the soliton defect is enhanced with the increase of U from 0.0 eV to 0.9 eV 
and then reduced substantially with the increase of U from 0.9 eV to 4.0 eV. This sequence 
is in good accordance with the sequence for the stationary velocity of a charged soliton Vg 
illustrated in Fig. P 
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FIG. 14: The staggered bond order parameter r„ of a static charged soHton for several different V 
values, (ep = 5.0 x 10"^ U=0.0 eV) 

Secondly, we also study the influence of nearest neighbor electron-electron interactions 
V on the charged soliton transport process. The dependence of the stationary velocity of 
a charged soliton Vs on the V values is displayed in Fig. [121 where U values are supposed 
to be frozen to be zero. This assumption is not realistic because the on-site Coulomb 
interactions U are normally much stronger than the neighbor electron-electron interactions 
V; we make this assumption only for the purpose of studying of the influence of V on the 
soliton transport process without the effect of U. As can be seen in Fig. [T21 the situation 
is completely different from the case of an on-site Coulomb interaction: Vg goes through a 
minimum for flnite V, not a maximum. Actually, similar to the the case discussed above 
with U varying, the increase or decrease of Vg with V is also strongly related to the change 
of the delocalization abilities of the charged soliton defect. In Fig. [T3] we show the net 
charge distribution pictures of a static charged soliton calculated with different V values. 
It can be clearly seen that the nearest neighbor electron-electron interactions V will induce 
negative charge densities at the neighboring sites of the positive charge density and the 
charge polarization will become larger and larger with increasing V. For small V, larger 
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positive charge densities will be induced in the central part of a charged soliton defect while 
the induced negative charge densities at nearest neighbor sites are very small. Therefore 
the charged soliton defect tends to be more localized and consequently the charged soliton 
transport becomes more difficult. So, the velocity of the charged soliton Vg decreases while 
the value of V increases from 0.0 eV to 1.1 eV as displayed in Fig. [121 However, when 
nearest neighbor electron-electron interactions V increase further and become dominant, 
large charge polarization will be induced. Therefore, the electron-hole attractions between 
opposite charge densities at nearest-neighbor sites will contribute much more significantly 
to the soliton system and favor the hoppings of the accumulated electrons (or holes) in the 
central part of a charged soliton defect to the neighbor sites. This change leads to a more 
delocalized soliton defect and accordingly the increase of the velocity of the charged soliton 
Vg. In order to directly view the change of delocalization level of the charged soliton defect 
through geometrical pictures, we also show the staggered bond order parameter r„ of a static 
charged soliton calculated with different V values in Fig. [Mj. We make hyperbolic tangent 
function fittings for these curves and obtain the soliton width L for different V values: 7.95 
(V^ = 0.0 eV), 6.08 {V = 1.1 eV) and 6.19 {V = 1.5 eV). The change of the soliton width 
L shows the change of the delocalization of the charged soliton. This sequence is in good 
agreement with the sequence of the stationary velocity of a charged soliton Vs calculated for 
different V values illustrated in Fig. [T21 

Furthermore, we consider the realistic case of trans-PA in which both the on-site Coulomb 
interactions U and the nearest-neighbor interactions V are taken into account. Fig. [T5l shows 
the stationary velocity of a charged soliton Vs as a function of U for fixed values of ^ = 0.2 
eV, V = 0.5 eV and V = 1.0 eV. Generally, the Vs behavior with fixed V values is quite 
similar to that has been illustrated in Fig. with fixed V=0.0 eV. It can be found that Vs can 
be remarkably improved with increasing U values at first and then gradually decreases with 
U after a shallow maximum. This result can also be easily understood based on the simple 
analysis in HM. For a fixed value of V, as the on-site repulsions U increase and gradually 
overcome the nearest-neighbor interactions V, the electron (or hole) will hop more easily to 
the next site, so the stationary velocity of a charged soliton increases. On the other hand, 
when the on-site repulsions U increase further and become to be dominant, the lattice tends 
to be occupied by one electron per site. The soliton motion will be greatly prevented by U, 
so the stationary velocity of a charged soliton decreases. Additionally, the U value where Vs 
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FIG. 15: The stationary velocity of a charged sohton as a function of U for fixed values of V. 
{At = 0.05 fs, Cp = 5.0 X 10-^ Eq = 3.0 mV/A) 

shows a shallow maximum increases while V is increased. 

IV. SUMMARY AND CONCLUSION 

For a model conjugated polymer chain initially holding a charged soliton defect, described 
using the combined SSH-EHM model extended with an additional part for the influence of 
an external electric field, we have solved the time- dependent Schrodinger equation for the 
TT electrons and the equations of motion for the monomer displacements in a nonadiabatic 
way by virtue of the combination of the adaptive TDDMRG and classical MD. 

We find that after a smooth turn-on of the external electric field the charged soliton 
is accelerated at first and then moves with a constant velocity as one entity consisting of 
both the charge and the lattice deformation. During the entire time evolution process the 
geometrical distortion curve and charge distribution shape for the charged soliton defect 
always stay coupled and show no dispersion, implying that the soliton defect is an inherent 
feature and the fundamental charge carrier in trans-PA. The dependence of the stationary 
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velocity of a charged soliton Vs on the external electric field strength is also studied, and an 
ohmic region where increases linearly with the field strength is found. This ohmic region 
extends approximately from 6 mV/A to 12 mV/A for the case of U=2.0 eV and ^=1.0 
eV. Below 6 mV/A the soliton velocity increases nonlinearly and no saturation of Vg with 
increasing electric field strength is found. 

The infiuence of electron-electron interactions (both the on-site Coulomb interactions 
U and the nearest-neighbor interactions V) on charged soliton motion are investigated in 
detail. In general, small U values are beneficial to the solitonic motion because they have the 
tendency to push the accumulated electrons (or holes) in the defect central part to hop more 
easily to the neighbor site, and large U values are unfavorable for the soliton motion because 
the lattice tends to be occupied by one electron per site due to the domination of on-site 
repulsions. Meanwhile, small V values are not beneficial to the solitonic motion because 
they will induce a more localized defect distribution, and due to the induced large charge 
polarization accompanied with large nearest-neighbor attractions large V values favor the 
soliton transport. When U and V are considered at the same time, Vg can be significantly 
enhanced at the beginning stage of increasing U and then decreases gradually with U after 
a shallow maximum. 

Because our calculations by the nonadiabatic TDDMRG/MD method give interesting 
charged soliton transport pictures and reveal the relationship between electron-electron in- 
teractions and the transport of charged solitons while all relevant electron-electron and 
electron-phonon interactions are nearly fully taken into account, we consider that this com- 
bination of methods provides a new nonadiabatic dynamic simulation method for exploring 
the physical properties of electrons correlation systems with relevant lattice dynamics. 
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